source('00_util_scripts/mod_seurat.R')
source('00_util_scripts/mod_bplot.R')
library(data.table)
library(tictoc)
library(furrr)
plan('multisession', workers = 8)

# chemotaxis uniprot keyword --------
chemotax <- query_uniprot_keyword('KW-0145')
cytokine <- query_uniprot_keyword('KW-0202')

chemokine <- chemotax |>
  filter(symbol %in% cytokine$symbol) |>
  write_csv('trpm7_mucosa/uniprot.chemokine.csv')

chemokine <- read_csv('trpm7_mucosa/uniprot.chemokine.csv')

# import ---------
h5.path <- list.files('trpm7_mucosa/data/', pattern = 'h5', full.names = T)

## 8 worker furrr only use 1/3 time of purrr::map
tic()
mex.list <- h5.path |>
  future_map(Read10X_h5, .progress = T)
toc()

mex.name <- h5.path |> str_extract('GSM.+GEX') |>
  str_remove_all('GSM\\d+_|_GEX')

mex.list[[1]][1:5,1:5]

mex.tidy <- mex.list |>
  map2(mex.name, add_name_field, .progress = T)

mex.tidy[[1]][1:5,1:5]

## primary -------------
prmry <- mex.name |> str_detect('X31|PR8', negate = T)

mex.prmry <- mex.tidy[prmry][[1]] |>
  RowMergeSparseMatrices(mex.tidy[prmry][-1])

colnames(mex.prmry) <- colnames(mex.prmry) |>
  str_remove('-1')

meta.prm <- read_delim('trpm7_mucosa/data/scp_primary_metadata.txt')

meta.prm[-1,] |>
  write_csv('trpm7_mucosa/primary.meta.csv')

meta.prm <- read_csv('trpm7_mucosa/primary.meta.csv')
# 156572 cells

barc.prm <- meta.prm$NAME |>
  str_replace('\\.', '_') |>
  str_replace('E','M')

meta.prm |>
  mutate(.cell = str_replace(NAME, '\\.', '_') |> str_replace('E','M'),
         tissue = str_extract(.cell, 'RM|OM|LNG'),
         time = str_extract(.cell, 'D\\d+|Naive'),
         barcode = str_extract(.cell, '[A-Z]+$'),
         .keep = 'none') |>
  dplyr::count(time)

barc.prm[1:5]

barc.prm.raw <- colnames(mex.prmry) |>
  str_replace('D2', 'D02') |>
  str_replace('D5', 'D05') |>
  str_replace('D8', 'D08')

colnames(mex.prmry) <- barc.prm.raw

barc.prm.int <- barc.prm.raw |>
  intersect(barc.prm)

sobj.prm <- mex.prmry[, barc.prm.int] |>
  CreateSeuratObject(min.cells = 3, min.features = 500)

sobj.prm |>
  write_rds('trpm7_mucosa/iav.prm.rds')

## recall ------
rech <- mex.name |> str_detect('X31|PR8|RM')

mex.re <- mex.tidy[rech][[1]] |>
  RowMergeSparseMatrices(mex.tidy[rech][-1])

sobj.re <- mex.re |>
  CreateSeuratObject(min.cells = 3, min.features = 500)

mex.re <- sobj.re |>
  GetAssayData()

mex.re[1:5,1:5]

colnames(mex.re) <- colnames(mex.re) |>
  str_remove('-1')

meta.re <-
  read_delim('trpm7_mucosa/data/scp_rechallenge_combined_metadata.txt')

meta.re[-1,] |>
  write_csv('trpm7_mucosa/recall.meta.csv')

meta.re <- read_csv('trpm7_mucosa/recall.meta.csv')
# 143330 cells

barc.re <- meta.re$NAME |>
  str_replace('\\.', '_') |>
  str_replace('E','M')

meta.re <- meta.re |>
  mutate(.cell = str_replace(NAME, '\\.', '_') |> str_replace('E','M'),
         barc = str_extract(.cell, '[A-Z]+$'),
         time = str_extract(.cell, 'D\\d+|Naive|C\\d'),
         group = case_when(str_detect(.cell, 'Naive') ~ 'Naive',
                           str_detect(.cell, 'RM') ~ 'Primary',
                           .default = str_extract(.cell, 'X31|PR8')),
         .cell = str_c(group, '.', time, '_', barc))

rework <- tibble(.cell = colnames(mex.re)) |>
  mutate(barc = str_extract(.cell, '[A-Z]+$'),
         time = str_extract(.cell, 'D\\d+|Naive|C\\d') |>
           str_replace('D2', 'D02') |>
           str_replace('D5', 'D05') |>
           str_replace('D8', 'D08'),
         group = case_when(str_detect(.cell, 'Naive') ~ 'Naive',
                           str_detect(.cell, 'RM') ~ 'Primary',
                           .default = str_extract(.cell, 'X31|PR8')),
         .cell = str_c(group, '.', time, '_', barc))

rework |>
  distinct(time, group)

barc.int <- rework$.cell |>
  intersect(meta.re$.cell)

rework$.cell[1:5]

mex.re |> colnames() |> head()

colnames(mex.re) <- rework$.cell

sobj.re <- mex.re[,barc.int] |>
  CreateSeuratObject(min.cells = 3, min.features = 500)

sobj.re |>
  write_rds('trpm7_mucosa/iav.rm.re.rds')

# rds preprocess ---------
## recall -----------
sobj.re <- read_rds('trpm7_mucosa/iav.rm.re.rds')

sobj.re$mito.ratio <- sobj.re |>
  PercentageFeatureSet('^mt-')

sobj.re |> VlnPlot('mito.ratio', pt.size = 0)

sobj.re <- sobj.re |>
  NormalizeData()

meta.re <- meta.re |>
  mutate(.cell, NAME, 
         label.main = cell_type__ontology_label,
         label.fine = cell_type_custom, .keep = 'none')

umap.re <- read_csv('trpm7_mucosa/umap.recall.csv')

# avoid using tidyseurat keywords
meta.re <- umap.re |>
  mutate(NAME, umap_1.old = X, umap_2.old = Y,
         label.coarse = Cell.Type, .keep = 'none') |>
  right_join(meta.re)

# 10 Cell.Type, 31 ctol, 126 custom label
sobj.re <- sobj.re |>
  left_join(meta.re)

sobj.re |>
  ggplot(aes(umap_1, umap_2, color = label.coarse)) +
  geom_point(size = .3) +
  scale_color_brewer(palette = 'Paired') +
  theme_classic() +
  labs(title = 'Cell types in respiratory mucosa')

sobj.re <- sobj.re |>
  mutate(time = fct_relevel(time, 'Naive') |>
           fct_relevel('C2','C5', after = Inf),
         orig.ident = fct_reorder(orig.ident, as.integer(time)) |>
           fct_relevel('X31.D60','X31.C2','X31.C5', after = Inf))

time.fct <- sobj.re |>
  distinct(time, orig.ident) |>
  mutate(time.chr = as.character(time),
         sample.chr = as.character(orig.ident), .keep = 'used')

## keep only epi + immune ---------
sobj.ei <- sobj.re |>
  filter(str_detect(label.coarse, 'B|Epi|Gran|Mye|T'))

sobj.ei |>
  ggplot(aes(umap_1, umap_2, color = label.coarse)) +
  geom_point(size = .5) +
  theme_classic() +
  labs(title = 'Cell types in respiratory mucosa')

sobj.ei <- sobj.ei |>
  quick_process_seurat(skip_norm = T)

new.umap <- sobj.ei |> Reductions('umap') |> pluck('cell.embeddings')

new.umap$umap_1 |> head()

sobj.ei <- sobj.ei |>
  mutate(umap_1.new = new.umap[,1], umap_2.new = new.umap[,2])

sobj.ei |> DimPlot(group.by = 'label.coarse') + ggtitle('')

### naive cell proportion --------
prop.ei <- sobj.ei |>
  filter(str_detect(orig.ident, 'Nai|D02')) |>
  calc_frac_conf_on_grouped_count(orig.ident, label.coarse)

prop.ei |>
  mutate(label.coarse = fct_relevel(label.coarse, 'Epithelial', 'Myeloid',
                                    'T/NK Cell', 'Granulocyte'),
         orig.ident = str_remove(orig.ident, 'Naive.')) |>
  ggplot(aes(label.coarse, fraction*100, fill = orig.ident)) +
  geom_col(position = 'dodge2') +
  theme_pubr(legend = 'right') +
  scale_fill_hue(direction = -1) +
  labs(title = 'Cell proportion in nasal respiratory mucosa',
       y = 'Cell proportion (%)', x = 'Cell type', fill = 'Group')

### expr in coarse type ----------
sobj.ei |>
  filter(str_detect(orig.ident, 'Nai')) |>
  mutate(`Cell type` = fct_relevel(label.coarse, 'Epithelial', 'T/NK Cell',
                                   'Granulocyte', 'Myeloid')) |>
  bill.violin('Trpm7', `Cell type`) +
  NoLegend() +
  labs(title = 'Trpm7 expression in naive nasal respiratory mucosa',
       y = 'Normalized expression')

m7.ei <- last_plot() |> pluck('data')

m7.ei |> write_csv('trpm7_mucosa/results/mice.naive.rm.epi.immu.m7.expr.csv')

m7.ei |>
  summarise(logtpm.mean(.abundance_RNA), .by = label.coarse)

m7.crs <- sobj.re |>
  DotPlot2d('Trpm7', label.coarse, group.y = orig.ident) |>
  pluck('data')

m7.crs |>
  filter(str_detect(group.y, 'Nai|Pri'), group.x != 'HSC') |>
  mutate(sample.chr = group.y,
         group.x = fct_reorder(group.x, avg.exp, .desc = T)) |>
  left_join(time.fct) |>
  ggplot(aes(group.x, orig.ident, size = pct.exp, color = avg.exp)) +
  geom_point() +
  scale_color_distiller(palette = "RdYlBu") +
  theme_pubr(legend = "right") +
  RotatedAxis() +
  labs(x = 'Cell type', y = 'Sample',
       color = "Average expression", size = "Percent expressed",
       title = 'Trpm7: RM primary IAV infection')

m7.crs |>
  filter(!str_detect(group.y, 'Nai|Pri'), group.x != 'HSC') |>
  mutate(sample.chr = group.y,
         group.x = fct_reorder(group.x, avg.exp, .desc = T)) |>
  left_join(time.fct) |>
  ggplot(aes(group.x, orig.ident, size = pct.exp, color = avg.exp)) +
  geom_point() +
  scale_color_distiller(palette = "RdYlBu") +
  theme_pubr(legend = "right") +
  RotatedAxis() +
  labs(x = 'Cell type', y = 'Sample',
       color = "Average expression", size = "Percent expressed",
       title = 'Trpm7: RM IAV rechallenge')

### expr in finer type --------
m7.main <- sobj.re |>
  DotPlot2d('Trpm7', label.main, group.y = orig.ident) |>
  pluck('data')

m7.main |>
  filter(str_detect(group.y, 'Nai|Pri')) |>
  mutate(sample.chr = group.y,
         group.x = fct_reorder(group.x, avg.exp, .desc = T)) |>
  left_join(time.fct) |>
  ggplot(aes(group.x, orig.ident, size = pct.exp, color = avg.exp.scaled)) +
  geom_point() +
  scale_color_distiller(palette = "RdYlBu") +
  theme_pubr(legend = "right") +
  RotatedAxis() +
  labs(x = 'Cell type', y = 'Sample',
       color = "Average expression", size = "Percent expressed",
       title = 'Trpm7: RM primary IAV infection')

m7.main |>
  filter(!str_detect(group.y, 'Nai|Pri')) |>
  mutate(sample.chr = group.y,
         group.x = fct_reorder(group.x, avg.exp, .desc = T)) |>
  left_join(time.fct) |>
  ggplot(aes(group.x, orig.ident, size = pct.exp, color = avg.exp.scaled)) +
  geom_point() +
  scale_color_distiller(palette = "RdYlBu") +
  theme_pubr(legend = "right") +
  RotatedAxis() +
  labs(x = 'Cell type', y = 'Sample',
       color = "Average expression", size = "Percent expressed",
       title = 'Trpm7: RM IAV rechallenge')

## epithelial -----------
sobj.epi <- sobj.re |>
  filter(label.coarse == 'Epithelial')

sobj.epi.rm.pr <- sobj.epi |>
  filter(str_detect(orig.ident, 'Naive|Pri'))

pri.time <- sobj.epi.rm.pr$orig.ident |>
  unique() |> str_subset('Pri')

pri.time |>
  map(\(x)FindMarkers(sobj.epi.rm.pr, features = 'Trpm7',
                      group.by = 'orig.ident', ident.1 = x,
                      ident.2 = 'Naive.Naive'))

sobj.epi.rm.pr |>
  ggplot(aes(umap_1, umap_2, color = label.main)) +
  geom_point(size = .5) +
  scale_color_brewer(palette = 'Paired') +
  theme_classic() +
  labs(title = 'Epithelial cells in respiratory mucosa')

sobj.epi.rm.pr |>
  DotPlot2d('Trpm7', label.main, orig.ident) +
  RotatedAxis()

sobj.epi.rm.pr |>
  filter(orig.ident == 'Naive.Naive') |>
  bill.violin('Trpm7', label.main) +
  NoLegend() + RotatedAxis()

epi.naive.m7 <- last_plot() |>
  pluck('data') 

epi.naive.m7 |>
  summarise(avg_expr = logtpm.mean(.abundance_RNA), .by = label.main)

Idents(sobj.epi.rm.pr) <- 'label.fine'

sobj.epi.rm.pr |>
  filter(orig.ident == "Naive.Naive") |>
  FindAllMarkers(features = 'Trpm7')

sobj.erp <- sobj.epi.rm.pr |>
  select(-c(umap_1, umap_2)) |>
  quick_process_seurat(skip_norm = T)

sobj.erp |>
  DotPlot2d('Trpm7', seurat_clusters, orig.ident) +
  RotatedAxis()

sobj.erp |>
  filter(orig.ident == "Naive.Naive") |>
  FindAllMarkers(features = 'Trpm7')

sobj.erp |>
  filter(seurat_clusters == 4) |>
  dplyr::count(label.fine, sort = T)

sobj.erp |>
  DimPlot(group.by = 'label.main', cols = 'Paired') +
  labs(title = 'Epithelial cells in respiratory mucosa')

### naive frac & trpm2-expr -----------
sobj.d02.epi <- sobj.re |>
  filter(str_detect(time, 'Naive|D02'), label.coarse == 'Epithelial')

epith_count <- sobj.d02.epi |>
  dplyr::count(label.main, time)

epith_count |>
  mutate(label.main = fct_reorder(label.main, n)) |>
  ggplot(aes(n, label.main, fill = label.main)) +
  geom_col() +
  facet_grid(~time) +
  labs(y = 'Cell type', title = 'Epithelial cell number in respiratory mucosa') +
  theme_bw()

epith_fine_count <- sobj.d02.epi |>
  dplyr::count(label.fine, time)

epith_fine_count |>
  mutate(label.fine = fct_reorder(label.fine, n)) |>
  ggplot(aes(n, label.fine, fill = label.fine)) +
  geom_col() +
  facet_grid(~time) +
  labs(y = 'Cell type', title = 'Epithelial cell number in respiratory mucosa') +
  theme_bw() + NoLegend()

sobj.d02.epi |>
  DotPlot2d('Trpm7', time, label.main) |>
  pluck('data') |>
  mutate(group.x = fct_relevel(group.x, 'Naive'),
         group.y = fct_reorder(group.y, avg.exp, max)) |>
  BubblePlot(d2 = T, size = c(1,10)) +
  labs(title = 'Trpm7 expression in RM epithelial cells',
       x = 'Time', y = 'Cell type')

sobj.d02.epi |>
  DotPlot2d('Trpm7', time, label.fine) |>
  pluck('data') |>
  mutate(group.x = fct_relevel(group.x, 'Naive'),
         group.y = fct_reorder(group.y, avg.exp, max)) |>
  BubblePlot(d2 = T, size = c(1,6)) +
  labs(title = 'Trpm7 expression in RM epithelial cells',
       x = 'Time', y = 'Cell type')

### rds checkpoint ---------
sobj.erp |>
  write_rds('trpm7_mucosa/rm.primary.epi.rds')

sobj.erp <- read_rds('trpm7_mucosa/rm.primary.epi.rds')

### pri vs naive dega ----------
epi.pri.fc <- sobj.erp |>
  FindMarkersAcrossVar(split.by = 'label.main', group.by = 'orig.ident',
                       ident.2 = 'Naive.Naive', ident.1 = 'Primary.D02')

epi.pri.fc |>
  write_csv('trpm7_mucosa/rm.primary.epi.d2vsNaive.deg.csv')

imm.pri.fc <- sobj.ei |>
  FindMarkersAcrossVar(split.by = 'label.coarse', group.by = 'orig.ident',
                       ident.2 = 'Naive.Naive', ident.1 = 'Primary.D02')

imm.pri.fc |>
  write_csv('trpm7_mucosa/epi.immune.d2vsNaive.deg.csv')

epi.pri.fc |>
  filter(gene %in% c('S100a8', 'S100a9')) |>
  ggplot(aes(cluster, gene, size = -log10(p_val_adj), fill = avg_log2FC)) +
  geom_point(shape=21) +
  scale_fill_gradient2(low = 'blue', high = 'red', limits = c(NA, 1)) +
  theme_bw() +
  RotatedAxis() +
  labs_pubr() +
  labs(x = 'Cell type')

imm.pri.fc |>
  filter(gene %in% str_to_title(chemokine$symbol)) |>
  mutate(gene = fct_reorder(gene, p_val_adj, .fun = min, .desc = T),
         cluster = fct_reorder(cluster, p_val_adj, .fun = min)) |>
  ggplot(aes(cluster, gene, size = -log10(p_val_adj), fill = avg_log2FC)) +
  geom_point(shape=21) +
  scale_fill_gradient2(low = 'blue', high = 'red') +
  theme_bw() +
  RotatedAxis() +
  labs_pubr() +
  labs(x = 'Cell type',
       title = 'Chemokine fold change 2 dpi vs Naive')

epi.pri.fc |>
  filter(gene %in% str_to_title(chemokine$symbol)) |>
  mutate(gene = fct_reorder(gene, p_val_adj, .fun = min, .desc = T),
         cluster = fct_reorder(cluster, p_val_adj, .fun = min)) |>
  ggplot(aes(cluster, gene, size = -log10(p_val_adj), fill = avg_log2FC)) +
  geom_point(shape=21) +
  scale_fill_gradient2(low = 'blue', high = 'red') +
  theme_bw() +
  RotatedAxis() +
  labs_pubr() +
  labs(x = 'Cell type',
       title = 'Chemokine fold change 2 dpi vs Naive:\nepithelial cells')

basal.gsego <- epi.pri.fc |>
  filter(cluster == 'basal cell', p_val_adj < .05) |>
  pull(avg_log2FC, name = gene) |>
  sort(decreasing = T) |>
  gseGO(OrgDb = 'org.Mm.eg.db', keyType = 'SYMBOL', eps = 0, ont = 'ALL')

basal.gsego <- basal.gsego |> simplify()

basal.gsego@result |>
  filter(!str_detect(Description, 'synapse')) |>
  plot_enrichment(metric = NES, n = 10) +
  labs(title = 'GSEA of RM basal epithelial cells after primary infection',
       x = 'NES')

## primary ---------
sobj.pr <- read_rds('trpm7_mucosa/iav.prm.rds')

umap.pr <- read_csv('trpm7_mucosa/primary.umap.csv')

umap.pr |> distinct(Cell.Type)

meta.pr <- meta.pr |>
  mutate(NAME, 
         label.main = cell_type__ontology_label,
         label.fine = cell_type_custom, .keep = 'none')

meta.pr <- umap.pr |>
  mutate(NAME, umap_1 = X, umap_2 = Y, label.coarse = Cell.Type, .keep = 'none') |>
  right_join(meta.pr)

meta.pr |>
  distinct(label.coarse, label.main) |>
  DT::datatable()

crs.mn <- meta.re |>
  dplyr::count(label.coarse, label.main) |>
  filter(n > 10)

meta.pr <- meta.pr |>
  select(-label.coarse) |>
  left_join(crs.mn) |>
  mutate(label.coarse = case_when(
    label.main == 'plasma cell' ~ 'B Cell',
    str_detect(label.main, '^glandul') ~ 'Epithelial',
    .default = label.coarse
  ))

meta.pr <- meta.pr |>
  select(-n) |>
  mutate(.cell = str_replace(NAME, '\\.', '_') |>
           str_replace('E', 'M'))

sobj.pr <- sobj.pr |> NormalizeData()

sobj.pr <- sobj.pr |>
  left_join(meta.pr)

sobj.pr |> write_rds('trpm7_mucosa/iav.prm.rds')

### OM -------------
sobj.om <- sobj.pr |>
  filter(str_detect(.cell, 'OM'))

sobj.om |>
  ggplot(aes(umap_1, umap_2, color = label.coarse)) +
  geom_point(size = .3) +
  scale_color_brewer(palette = 'Paired') +
  theme_classic() +
  labs(title = 'Cell types in olfactory mucosa')

#### expr in coarse type --------
m7.crs.om <- sobj.om |>
  DotPlot2d('Trpm7', label.coarse, group.y = orig.ident) |>
  pluck('data')

m7.crs.om |>
  mutate(group.y = fct_relevel(group.y, 'Naive'),
         group.x = fct_reorder(group.x, avg.exp, .desc = T)) |>
  ggplot(aes(group.x, group.y, size = pct.exp, color = avg.exp.scaled)) +
  geom_point() +
  scale_color_distiller(palette = "RdYlBu") +
  theme_pubr(legend = "right") +
  RotatedAxis() +
  labs(x = 'Cell type', y = 'Time',
       color = "Average expression", size = "Percent expressed",
       title = 'Trpm7: OM primary IAV infection')

#### expr in main type --------
m7.main.om <- sobj.om |>
  DotPlot2d('Trpm7', label.main, group.y = orig.ident) |>
  pluck('data')

m7.main.om |>
  mutate(group.y = fct_relevel(group.y, 'Naive'),
         group.x = fct_reorder(group.x, avg.exp, .desc = T)) |>
  ggplot(aes(group.x, group.y, size = pct.exp, color = avg.exp.scaled)) +
  geom_point() +
  scale_color_distiller(palette = "RdYlBu") +
  theme_pubr(legend = "right") +
  rotate_x_text(45) +
  labs(x = 'Cell type', y = 'Time',
       color = "Average expression", size = "Percent expressed",
       title = 'Trpm7: OM primary IAV infection')

#### expr in fine type --------
m7.fine.om <- sobj.om |>
  DotPlot2d('Trpm7', label.fine, group.y = orig.ident) |>
  pluck('data')

m7.fine.om |>
  mutate(group.y = fct_relevel(group.y, 'Naive'),
         group.x = fct_reorder(group.x, avg.exp, .desc = T)) |>
  slice_min(group.x, n = 150) |>
  ggplot(aes(group.x, group.y, size = pct.exp, color = avg.exp)) +
  geom_point() +
  scale_color_distiller(palette = "RdYlBu") +
  theme_pubr(legend = "right") +
  rotate_x_text(angle = 45) +
  labs(x = 'Cell type', y = 'Time',
       color = "Average expression", size = "Percent expressed",
       title = 'Trpm7: OM primary IAV infection')

### LNG ------------
sobj.pr |>
  filter(str_detect(.cell, 'LNG')) |>
  ggplot(aes(umap_1, umap_2, color = label.coarse)) +
  geom_point(size = .3) +
  scale_color_brewer(palette = 'Paired') +
  theme_classic() +
  labs(title = 'Cell types in lateral nasal gland')

sobj.lng <- sobj.pr |>
  filter(str_detect(.cell, 'LNG'))

#### expr in coarse type --------
m7.crs.lng <- sobj.lng |>
  DotPlot2d('Trpm7', label.coarse, group.y = orig.ident) |>
  pluck('data')

m7.crs.lng |>
  mutate(group.y = fct_relevel(group.y, 'Naive'),
         group.x = fct_reorder(group.x, avg.exp, .desc = T)) |>
  ggplot(aes(group.x, group.y, size = pct.exp, color = avg.exp.scaled)) +
  geom_point() +
  scale_color_distiller(palette = "RdYlBu") +
  theme_pubr(legend = "right") +
  RotatedAxis() +
  labs(x = 'Cell type', y = 'Time',
       color = "Average expression", size = "Percent expressed",
       title = 'Trpm7: LNG primary IAV infection')

#### expr in main type --------
m7.main.lng <- sobj.lng |>
  DotPlot2d('Trpm7', label.main, group.y = orig.ident) |>
  pluck('data')

m7.main.lng |>
  mutate(group.y = fct_relevel(group.y, 'Naive'),
         group.x = fct_reorder(group.x, avg.exp, .desc = T)) |>
  ggplot(aes(group.x, group.y, size = pct.exp, color = avg.exp.scaled)) +
  geom_point() +
  scale_color_distiller(palette = "RdYlBu") +
  theme_pubr(legend = "right") +
  rotate_x_text(60) +
  labs(x = 'Cell type', y = 'Time',
       color = "Average expression", size = "Percent expressed",
       title = 'Trpm7: LNG primary IAV infection')
